#!/bin/csh
# 1d_profiles.csh
################################################################
# Plots of horizontally averaged profiles of different fields
# 
# - Set the simulation directory containing YACC's output 
#   in the file "simulation.dir"
# - Run the script using:
#   ./1d_profiles.csh n 
#   'n' is the time step to be plotted 
#    Note that n must be given in the following form:
#    000001 or 000010 or 000100 or 001000 etc.
###########################################################

#-----------------------------------------------------------------------------------
# Possible profiles:
#-----------------------------------------------------------------------------------
# temperature (default: 1)
# velocity: velocity yes or no (default: 1)
# urms: RMS velocity (default: 1)
# ux: horizontal velocity (default: 1)
# uz: vertical velocity (default: 1)
# viscosity (default 0)
# viscdiss: viscous dissipation (default 0)
# composition (default: 0)
# heats : Rayleigh number for internal heating (default: 0)
# conductivity: thermal conductivity (default: 0)
# expansivity: thermal expansivity (default: 0)
# solidus0: initial solidus (default: 0)
# solidus: current solidus (default: 0)
# liquidus (default: 0)
# depletion (default: 0)
# surfvel: horizontal profile of the surface velocity (default: 0)
#
# Select profiles to plot
set temperature = 1
set velocity = 1
set urms = 1
set ux = 1
set uz = 1
set viscosity = 1
set viscdiss = 0
set composition = 0
set heats = 0
set conductivity = 0
set expansivity = 0
set depletion = 0
set solidus0 = 0
set solidus = 0
set liquidus = 0
set surfvel = 0
#-----------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------
# Options:
#-----------------------------------------------------------------------------------
# dirout: directory of the YACC's output files (needs to be written in the 
#         'simulation.dir' file)
# dirplots: directory where to store the plots (default: $dirout/Plots/$field)
# makejpg: make also jpg files (default: 1)
# jpgdensity: resolution of the converted jpg files (default: 150)
# deletepdf: delete pdf files (done only if makejpg=1, default: 1)
# dimx: plot width in cm (default: 10 cm)
# dimz: plot height in cm (default: 10 cm)
# X0: x-coordinate in cm of the bottom-left corner of the plot (default: 3 cm)
# Y0: Y-coordinate in cm of the bottom-left corner of the plot (default: 5 cm)
# zmin: z-coordinate of the bottom edge of the plot (default: 0)
# zmax: z-coordinate of the top edge of the plot (default: 1)
#
# Select options
set makejpg = 1
set jpgdensity = 150
set deletepdf = 1
set dimx = 10
set dimz = 10
set dimx2 = 15
set dimz2 = 5
set X0 = 3
set Y0 = 5
set zmin = 0
set zmax = 1
#-----------------------------------------------------------------------------------

#--------------------------------------------------------------------------------------
# Specific options are set separately for each field (see below):
#--------------------------------------------------------------------------------------
# xlabel: label for x-axis (default: field) 
# zlabel: label for z-axis (default: Depth)
# min: set manually the minimum of the field 
# max: set manually the maximum of the field
# deltax: annotation interval for x-axis (default: (max-min)/4)
# deltaz: annotation interval for z-axis (default: 0.2)
# pdfoutput: name of pdf output file (default field_profile_{$tstep}.pdf)
# jpgoutput: name of jpg output file (default field_profile_{$tstep}.jpg)
#--------------------------------------------------------------------------------------

source simulation.dir

if (! -d $dirout) then
  echo $dirout 'does not exist. Check the path stored in the variable "dirout" '
  exit
endif

#------------------------------
# Set ouput directories
#------------------------------
set dirplots0 = {$dirout}Plots/
set dirplots = {$dirout}Plots/Profiles/ 

#-------------------------------
# Select YACC's profile file
#-------------------------------
if (! -d $dirplots0) then
  mkdir $dirplots0
endif

if (! -d $dirplots) then
  mkdir $dirplots
endif

set tstep = ${1} 
set nchar = `echo $tstep | awk '{print length($0)}'`
if ($nchar != 6) then
  echo "Wrong format of the input time step"
  exit
endif

set gridfile = ${dirout}grd_{$tstep}_.dat
set profiles = ${dirout}profiles_{$tstep}_.dat
echo "YACC grid output file:" $gridfile
echo "YACC profile output file:" $profiles

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

#------------------------------------------------------------
# Read and set current simulation time in the variable ctime
# (it can be used e.g. as plot title)
#------------------------------------------------------------
set ctime = `awk 'NR == 1 {print $5}' $gridfile`

####################
# Temperature
####################
if ($temperature == 1) then
  echo "Temperature profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Temperature"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = temperature_profile_{$tstep}.pdf
  set jpgoutput = temperature_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$3}; if($3< min) {min=$3}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$3}; if($3>max) {max=$3}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $3, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Velocity
####################
if ($velocity == 1) then
  echo "Velocity profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Velocity"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = velocity_profile_{$tstep}.pdf
  set jpgoutput = velocity_profile_{$tstep}.jpg

  set min = 0
  set max = `awk 'NR > 1 {if(min==""){min=max=$2}; if($2>max) {max=$2}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=0; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  # RMS 
  if ($urms == 1) then
    awk 'NR>1 {print $2, 1-$1}' $profiles| \
    psxy -R -JX -W7/0/0/0 -O -K >> $psoutput
  endif
  # Horizontal 
  if ($ux == 1) then
    awk 'NR>1 {print $7, 1-$1}' $profiles| \
    psxy -R -JX -W7,0/0/255 -O -K >> $psoutput
  endif
  # Vertical 
  if ($uz == 1) then
    awk 'NR>1 {print $8, 1-$1}' $profiles| \
    psxy -R -JX -W7,255/0/0 -O >> $psoutput
  endif

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Surface velocity
####################
if ($surfvel == 1) then
  echo "Surface velocity for time step" $tstep "at time" $ctime

# Read aspect ratio
  set ar = `awk 'NR == 1 {print $11}' $gridfile`
  if ($ar == "") then
    echo 'You are using an old output file with no information on the aspect ratio.'
    echo 'Set the aspect ratio explicitly in 2d_plots.csh as follows: "set ar = ... "'
    exit
  endif

# Read number of vertical points
  set nz = `awk 'NR == 1 {print $9}' $gridfile`
  set nz1 = `echo "scale=0; $nz-1" | bc -l`

  set title = ""
  set xlabel = "Horizontal@~\040@~distance"
  set zlabel = "Surface@~\040@~velocty"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = surfvelocity_profile_{$tstep}.pdf
  set jpgoutput = surfvelocity_profile_{$tstep}.jpg

# Store the x-coordinate and a the surface horizontal velocity in temporary file p_$tstep
  awk -v nz=`echo $nz` -v nz1=`echo $nz1` 'NR%nz==nz1 {print $1, $3}' $gridfile > p_$tstep

  set min = `awk '{if(min==""){min=$2}; if($2< min) {min=$2}} END {printf "%12.5f", min}' p_$tstep`
  set max = `awk '{if(min==""){min=max=$2}; if($2>max) {max=$2}} END {printf "%12.5f", max}' p_$tstep`
  set deltax = 1
  set deltaz = `echo "scale=1; ($max - $min)/4" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R0/$ar/$min/$max/ -JX$dimx2/$dimz2 -Ba{$deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk '{print $1, $2}' p_$tstep | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Viscosity
####################
if ($viscosity == 1) then
  echo "Viscosity profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "log@-10@-(Viscosity)"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = viscosity_profile_{$tstep}.pdf
  set jpgoutput = viscosity_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=log($4)/log(10)}; if(log($4)/log(10)< min) {min=log($4)/log(10)}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($4)/log(10)}; if(log($4)/log(10)>max) {max=log($4)/log(10)}} END {printf "%12.5f", max}' $profiles`
  if ($min == $max) then
    set min = -1
    set max = 1
  endif
  set deltax = `echo "scale=1; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $psoutput
  awk 'NR>1 {print log($4)/log(10), 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

######################
# Viscus dissipation
######################
if ($viscdiss == 1) then
  echo "Viscous dissipation profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "log@-10@-(Viscous@~\040@~dissipation)"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = viscous_dissipation_profile_{$tstep}.pdf
  set jpgoutput = viscdiss_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=log($5)/log(10)}; if(log($5)/log(10)< min) {min=log($5)/log(10)}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($5)/log(10)}; if(log($5)/log(10)>max) {max=log($5)/log(10)}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=1; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $psoutput
  awk 'NR>1 {print log($5)/log(10), 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

######################
# Composition
######################
if ($composition == 1) then
  echo "Composition profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Composition"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = composition_profile_{$tstep}.pdf
  set jpgoutput = composition_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$6}; if($6< min) {min=$6}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$6}; if($6>max) {max=$6}} END {printf "%12.5f", max}' $profiles`
  if ($min == $max) then
    set min = -1
    set max = 1
  endif
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $psoutput
  awk 'NR>1 {print $6, 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Conductivity
####################
if ($conductivity == 1) then
  echo "Conductivity profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Thermal@~\040@~conductivity"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = conductivity_profile_{$tstep}.pdf
  set jpgoutput = conductivity_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$11}; if($11< min) {min=$11}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$11}; if($11>max) {max=$11}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $11, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Expansivity
####################
if ($expansivity == 1) then
  echo "Expansivity profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Thermal@~\040@~expansivity"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = expansivity_profile_{$tstep}.pdf
  set jpgoutput = expansivity_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$10}; if($10< min) {min=$10}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$10}; if($10>max) {max=$10}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $10, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

######################
# Heat sources
######################
if ($heats == 1) then
  echo "Heat sources profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Ra@-H@-"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = heats_profile_{$tstep}.pdf
  set jpgoutput = heats_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$9}; if($9< min) {min=$9}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$9}; if($9>max) {max=$9}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=1; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $psoutput
  awk 'NR>1 {print $9, 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Initial solidus
####################
if ($solidus0 == 1) then
  echo "Initial solidus profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Inital@~\040@~solidus"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = solidus0_profile_{$tstep}.pdf
  set jpgoutput = solidus0_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$13}; if($13< min) {min=$13}} END {printf "%13.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$13}; if($13>max) {max=$13}} END {printf "%13.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $13, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Depletion
####################
if ($depletion == 1) then
  echo "Depletion profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Depletion"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = depletion_profile_{$tstep}.pdf
  set jpgoutput = depletion_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$12}; if($12< min) {min=$12}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$12}; if($12>max) {max=$12}} END {printf "%12.5f", max}' $profiles`
  if ($min == $max) then
    set min = -1
    set max = 1
  endif
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $12, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# Solidus
####################
if ($solidus == 1) then
  echo "Solidus profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Solidus"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = solidus_profile_{$tstep}.pdf
  set jpgoutput = solidus_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$14}; if($14< min) {min=$14}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$14}; if($14>max) {max=$14}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $14, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif

####################
# liquidus
####################
if ($liquidus == 1) then
  echo "Liquidus profile for time step" $tstep "at time" $ctime

  set title = ""
  set xlabel = "Liquidus"
  set zlabel = "Depth"

  set psoutput = temperature_profile_{$tstep}.ps
  set pdfoutput = liquidus_profile_{$tstep}.pdf
  set jpgoutput = liquidus_profile_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$15}; if($15< min) {min=$15}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$15}; if($15>max) {max=$15}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $psoutput
  awk 'NR>1 {print $15, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $psoutput

  ps2raster $psoutput -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

  if ($makejpg == 1) then
    mv -f $jpgoutput $dirplots
  endif
  if ($makejpg == 1 && $deletepdf == 1) then
    rm -f $pdfoutput
  else 
    mv -f $pdfoutput $dirplots
  endif
  rm -f p_$tstep $psoutput 

endif


